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Solutions to the random Fibonacci recurrence x n +i = x n ± j3x n -\ decrease (increase) exponen- 
tially, x n ~ exp(An), for sufficiently small (large) f3. In the limits j3 — > and (5 — > oo, we expand the 
Lyapunov exponent X(f3) in powers of /3 and /3 _1 , respectively. For the classical case of (3 = 1 we 
obtain exact non-perturbative results. In particular, an invariant measure associated with Ricatti 
variable r n — x n+ i/x n is shown to exhibit plateaux around all rational r's. 
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I. INTRODUCTION 

The Fibonacci numbers 1,1,2,3,5,8,13,... defined via 
F n _|_i = F n + F n -\ abound in nature For example, 

they govern the number of leaves, petals and seed grains 
in plants ||; they also give the number of ancestors of 
a drone. The Fibonacci recurrence is the simplest recur- 
rence in which each number depends on the previous two 
and this perhaps explains why the Fibonacci sequence is 
so ubiquitous. A natural stochastic modification of the 
Fibonacci sequence is to allow both additions and sub- 
stractions. Random Fibonacci sequences are related to 
many fields including condensed matter physics, dynam- 
ical systems, products of random matrices Q], continued 
fractions, etc. Random recurrences also form a chapter 
of the larger subject of iterated random functions ||. 

The random Fibonacci recurrence x n +\ = x n ± x n -\ 
results in sequences which behave erratically for small n. 
In the limit n — > oo, however, exponential growth occurs 
with probability one as was established by Furstenberg 
H in 1963. The large n behavior is characterized by the 
Lyapunov exponent A, 



A = lim 



ln\x n \ 



(1) 



Exponential growth (decay) means that A is positive 
(negative). For the random Fibonacci recurrence where 
each ± sign is independent and either + or — with prob- 
ability 1/2, the Lyapunov exponent is A = 0.12397559 . . . 
. Exponential growth may seem unsurprising but sim- 
ilar generalized random Fibonacci recurrence, 



X n _)_i — X n i pXyi—\) 



(2) 



gives exponential growth only when the parameter (3 
is sufficiently large, (3 > (3 S ~ 0.70258, whereas for 
< (3 < (3 S solutions decay exponentially Q. The de- 
cay occurs even though the expected values (x n ) remain 
constant and the expected values of the higher integer 
moments (a: 2 ), (x„), etc. grow exponentially. Indeed, 
from Eq. (||) one can deduce (xf l+1 ) = (x 2 ) + /3 2 (xf l _ 1 ) 
which implies (a; 2 ) + (j + /3 2 ) 1 / 2 ] , and similarly 

for higher moments M. Additionally, for (3 > 1/4 an 
interesting non-smooth dependence of the Lyapunov ex- 
ponent on the parameter (3 has been observed Q suggest- 
ing that the curve A(/3) is a fractal (it remains unclear 



whether this curve becomes genuinely smooth for suffi- 
ciently large (3). Similar non-smooth dependence on pa- 
rameters has been reported for several disordered systems 
0,0. Also, numerical results ||] suggest the following 
asymptotic behaviors of the Lyapunov exponent: 




when p 
when f3 



0, 



oo . 



(3) 



While understanding of the nature of the curve A(/3) 
might be a very difficult problem, one should be able to 
carry out asymptotic expansions of the Lyapunov expo- 
nent using tools developed in studies of one-dimensional 
disordered systems pp|-^7)|. Indeed, Eq. (||) admits the 
standard interpretation in terms of the one-dimensional 
(discretized) Schrodinger equation. A peculiarity of the 
present problem is that the corresponding random Hamil- 
tonian is non-Hermitian. Non-Hermitian random Hamil- 
tonians appear in various non-equilibrium problems |l§| l 
and exhibit interesting behaviors [|l9-22|, e.g., a der- 
ealization transition may occur even in one dimension. 
Fortunately, tools developed for Hermitian problems can 
often be extended to the non-Hermitian case. 

In the next Sect. II, we employ perturbation theory 
to expand the Lyapunov exponent in powers of (3. In 
particular, we show that in the [3 — > limit the second 
term is -1.75 x f3 4 rather than -1.875 x /3 4 [cf. Eq. (g)] 
and compute many more terms. In Sec. Ill, we use non- 
perturbative techniques to analyze the special f3 = 1 case. 
In Sec. IV, we examine random sequences with Gaussian, 
rather than binary, disorder. Finally, a few open ques- 
tions are discussed in the last Sec. V. 



II. PERTURBATION THEORY 



An old way of studying linear random recursions is to 
reduce them to random maps by introducing the Ricatti 
variable R n — x n+ i/x n . In the present case, Eq. (|^) 
becomes 



Rn = l± 



Rr, 



(4) 



1 



The Lyapunov exponent is given by 

1 " 

A = lim — } In Rj . 

n — >no n — ' 



(5) 



3=1 



In the large n limit, the probability distribution of R n 
becomes independent on n. The invariant distribution 
P{R,[3) satisfies 



P(R) 



2{R-1) 



P 



R-l 



P 



l-R 



(6) 



which immediately follows from the random map (Q). 
The analysis of the functional equation (JsJ) often sim- 
plifies after transforming it into an integral equation 



P(R,P) = J dR'P(R',/3) 



(7) 



8[R-l + jg 



Once we know the invariant distribution P{R, (3), we can 
employ Eq. (|^) to compute the Lyapunov exponent 



X(j3) = J dRP{R,f3) InR. 



(8) 



Depending on the magnitude of /3, the support of the 
invariant distribution P(R, (3) may be either finite or infi- 
nite. Let us assume the former. Then the extreme values 



satisfy R n 



R\x\\n — 



1 — fi/Rmin an d R n 



(3/Rn 



„ ; "max — „ ■ V y J 



Thus, the support is finite when the strength of disorder 
satisfies (3 < 1/4. Furthermore, in this case the sup- 
port is not merely the interval [i? m i n , ^?max] but rather 
a fractal set similar to the Cantor set. Indeed, the map 
R — ► 1 ± (3/R transforms the interval [i? m i n , i? ma x] into 
the union of two subintervals, [i? m i n , 1 — /3/-R max ] and 
[1 + /J/-R m ax,-Rmax]- Restricting the map to these two 
subintervals shows that they are transformed into four 
subintervals. Proceeding in the same manner ad infini- 
tum we construct the fractal support of the invariant dis- 
tribution. On the other hand, the support is unbounded 
when j3 > 1/4. This suggests to employ different pertur- 
bation approaches for small and large (3. 



A = 1 I dRP(R) 



In 1 - 



In 1 



(10) 



Expanding the logarithms on the right-hand side of 
Eq. ( |l(i| ) we obtain 

A = J dRR- 2 P{R) - J dRR- 4 P(R) + ... 

In the limit (3 — > 0, the interval [i? m in, -Rmax] shrinks to 
R = 1 . Hence, the first integral on the right-hand side of 
the above equation approaches to j dRP(R) = 1. Thus 
A = —f3 2 /2 + C(/3 4 ) and the first term of the expansion 
was indeed derived without using an explicit form of the 
invariant distribution. However, the derivation of the 
next term is still impossible without knowledge of the in- 
variant distribution. To avoid this we transform Eq. ( |Io| ) 
as we transformed Eq. (||) into Eq. ( |io| ) . Namely, we plug 
Eq. (@) into Eq. (0) to give 



A = 1 / ,1RP(R. /)!,(/?. /) 



(11) 



where 



L 2 (R,(3) =ln 1 



[3 



R , 



In 1 - 



1 + 1 



In 1 



l + £ 

' R , 



In 1 + 



1 - 



Expanding the logarithms and inserting this expansion 
into Eq. (|ll| ) we obtain 

A = -y / dRP(R)-^ J dR (± + l^P(R) + ... 

The first integral is now identically equal to one, while 
the second approaches to 7 in the small (3 limit. There- 
fore, A = -(3 2 /2 - 7/3 4 /4 + 0((3 6 ) is the two-term expan- 
sion. This shows that the small (3 expansion of Eq. (|3|) is 
slightly incorrect. 

Repeating the above procedure k times, we recast 
Eq. (|8|) into 



A 



dRP{R,/3) L k (R,0), 



(12) 



A. Weak disorder expansion 

When < 1/4, it is desirable to compute A(/3) with- 
out explicit determination of the invariant distribution 
P(R, (3) which is a complicated singular function. The 
trick is to transform the integral in the basic relation (||) 
into an integral which can be calculated perturbatively 
using only the normalization requirement J dR P(R) = 1. 
To this end we insert Eq. (0) into Eq. (||) to yield 



where 

L k (R,/3)= J2 ln[l;ei/3,l,e 2 /3,...,e fc /3,i?]. (13) 

The sum on the right-hand side of Eq. (|l^) is taken over 
all ej — ±1, and [1; a, 6, . . .] denotes the continued frac- 
tion 1 + 5^—- Expanding Lk{R,f3) and plugging it into 
Eq. ([l2]) one finds the correct expansion of the Lyapunov 
exponent up to 0((3 2k ). A few of these terms can be 



2 



computed by hand, and a bit more can be done with the 
help of Mathematics. One gets 



A(/3) = ~/? 2 -J/3 4 -f P 6 - 
2843 10 30755 12 
5 P 6 P 



555 



G(/3 14 ) 



(14) 



The radius of convergence of this series appears to be 
equal to 1/4 as one could guess from Eqs. (g). Hence the 
Lyapunov exponent is an analytic function of (3 when 
\(3\ < 1/4. Amusingly, the invariant distribution is a 
very singular function in this range of the parameter (3. 
The series (|l4|) perfectly reproduces numerical results ||] 
(5 representative digits), except for the case j3 = 1/4 for 
which A num w -0.0429 and A t i 1C or ~ -0.0424. This slight 
discrepancy is due to the fact that for (3 — 1/4 the generic 
term of the series decays algebraically in contrast with an 
exponential decay for (3 < 1/4. 



B. Strong disorder expansion 

In the large [3 limit, the support of the invariant dis- 
tribution P(R,f3) is the whole real line. The trick which 
has been employed in the case of weak disorder does not 
apply, i.e., we cannot determine A(/3) without knowledge 
of the invariant distribution. It proves convenient to use 
the modified Ricatti variable r n — \x n+ i/x n \fj3\. Then, 
Eq. (Q) reduces to the random map 



1 



±6 



8 = f3- 1 ' 2 . 



(15) 



Once we know the invariant distribution P(r, 5), we can 
compute the Lyapunov exponent via Eq. (|8|) which now 
becomes 



1 



A = - ln/3 + / drP(r) lnr. 



(16) 



The support of the invariant distribution is < r < oo. 
It proves convenient to define P(r) for negative r via 
P(r) = P(—r), so the support is the entire real line. In 
present notations, the functional equation (^) becomes 



2P(r) = 



(r + sy 



■P 



+ 5 J (r-5) 2 



P 



-S 



(17) 



For S — 0, this equation reduces to P(r) = r~ 2 P(l/r) 
which has infinitely many solutions. For S > 0, how- 
ever, the invariant distribution is unique. Thus taking 
the 5 — > limit of the invariant distribution P(r, S) we 
shall obtain a unique appropriate solution. In this sense, 
we are building a degenerate perturbation approach. 

To determine P(r, <5), notice that Eq. (|l7) can be re- 
written as 

P(r) = r- 2 P(l/r) + ~ D 2 [r- 2 P(l/r)] , (18) 



where D 2 is the discrete analog of the second derivative, 
D 2 F(r) = F(r + S) - 2P(r) + F{r - S). Expanding the 
right-hand side of Eq. ([18]) into Taylor series and denot- 
ing D = 4- we obtain 



P(r) = 



00 g2k 

^ (2fc)I 

k=0 v ; 



D 



2k \-2 



P(l/r)] 



(19) 



The change of variable r — * 1/r transforms Eq. ( p^ ) into 

(20) 



oo ?2k 



k=a 



where V = r 2 If we now plug Eq. @ into Eq. (|T|), 
we eliminate P(l/r) and thus arrive at a closed equation 
for the invariant distribution P(r, 5): 



oc oc 



J2fc+2Z 



D 2k { r -2 v 2l [ r 2p( r )] } . (21) 



Equation ( pT| ) suggests to seek a solution P{r, 8) which 
can be expanded as a power series in S 2 : 

oo 

P(r,5)=J2^P 3 (r). (22) 
j=o 

Plugging the series ( p2| ) into Eq. (^l|) one finds that the 
terms of order one cancel. Equating terms of order S 2 
leads to the following second order differential equation 
for Po(r) 



Dd P (r) = with £i = D + r 2 Dr 2 . 



(23) 



Integrating once we arrive at £iPo = 0, the integration 
constant being equal to zero to ensure that P (r) vanishes 
in the large r limit. This can be re- written as 



(l+r 4 )PoM + 2r 3 P (r) = 0, 



whose normalized solution is 



Po(r) 



20F 1 

r 2 (i/4) VTT^ 1 



(24) 



(25) 



Thus among infinitely many invariant distributions sat- 
isfying the duality relation P(r) = r~ 2 P(\/r) we have 
indeed selected the specific solution (p5|). Note that 
/ dr P (r)\nr = (this is actually valid for any func- 
tion which obeys the duality relation). Therefore, the 
anticipated contribution of order one to the Lyapunov 
exponent [cf. Eq. (0)] is actually equal to zero. 

Similarly equating the terms of order S 4 one finds 
DC\P\ +DC3P0 = 0, where £3 is certain third order dif- 
ferential operator. Integrating the above equation yields 
C\P\ + C3P0 = (the integration constant is equal to 
zero to ensure that P\{r) vanishes in the large r limit). 
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Recalling that CiPq = one can simplify the term C 3 P$. 
The final first order differential equation for Pi (r) reads 



(5 + r 4 )J 3 + 2r 3 g 2 
£>\Pi{r) = ^ Po(r). 



(26) 



It is helpful to write Pi (r) = Po (r) fi (r) and then ex- 
ploit the identity Ci[P (r) fi(r)] = (1 + r 4 )P {r)f[{r). 
Integrating the resulting equation yields 



fi(r)=ci+Ri(r), Ri(r) 



! (3r 8 + 8r 4 -15) 
6(r 4 + l) 3 ' 



(27) 



The integration constant ci can be found from the rela- 
tion J dr Pi (r) = which ensures that the normalization 
condition J dr P(r) — 1 holds. We get 



ci 



£(3/4) 

Lr(i/4) 



~] 2 



(28) 



A direct computation gives J dr Pi(r) hir = c\. Thus, in 
this order the contribution to the Lyapunov exponent is 
equal to ci/3 -1 . 

To determine Pj(r), we repeat the above procedure. 
Plugging into (|2l|), equating terms of order <5 2j+2 , 
and integrating once the resulting equation yields 



CiP,{r) + £ C 2n +iPj-n(r) = 0, 



(29) 



n=l 



where Cm+\ is the differential operator of order 2n + 1, 
2 



(2n + 2)! 

n 



2 



fc=0 



(2fc)!(2n - 2fc + 2) 



£) 2fc252(ri-fc) + l 7 ,2 



Equations (g^) can be solved recursively. Writing again 
Pj(r) — fj(r)P$(r), we find in the second order 



/ drPj(r) — and the basic identity T(l + z) = zT(z) 
for the gamma function |2|| . 

After inserting above expressions for Pj(r) in Eq. (|l|), 
we arrive at the expansion 



00 

A(/3) = -lnj8 + 5^a fcJ 9- 



(30) 



fc=i 



One can in principle generate A(/3) to any order. The 
first few coefficients are 



ax = ci = 0.1142366452611159 . . . 
0,2 = 4 + — = 0.01999445564958. . 



144 

a 3 = 0.0105345239 . 
a 4 = 0.0176632096. 



The coefficients 's have extremely complicated analyti- 
cal expressions in terms of the T function which we have 
not been able to simplify (although it appears plausi- 
ble that these expressions can be reduced to polynomials 
of ci with rational coefficients). The exact value for ai 
is in good agreement with the numerical estimate from 
Rcf. From the first four coefficients, one could guess 
that the radius of the convergence of series (|3(]) is of or- 
der unity. However, an apparent fractal structure of the 
curve A(/3) makes the above guess questionable. Even if 
our strong disorder expansion is asymptotic, it is quite 
accurate as can be seen by comparison of the four-term 
expansion (^0|) with numerical results. For (3 = 8, both 
give A « 1.05433 whereas there is a slight discrepancy 
for (3 = 4 as A num » 0.72309 and A th cor « 0.72319. 
On Fig. 1, we compare analytical (third-order expansion) 
and numerical results for the invariant distribution. 



h{r)=-- + cih{r)+R 2 {r), 



Mr) 



r 4 (463 - 3640r 4 + 2514r 8 + 440r 12 - 17r 



16 \ 



24(r 4 + l) 6 



In the third order we have 
1 



, , , 14699 



21600 96 

r 2 Q 3 (r) 



fi(r) + cif 2 (r)+R 3 (r) 



= 15120(r 4 + l) 9 
Q 3 (r) = 11340r 32 - 67 8 8 2 5r 28 - 112 6 3 68r 24 

- 3619377r 20 + 35 6 87 1 2 72r 16 - 471736467r 12 
+ 125696592r 8 - 5587155r 4 + 11340. 

The constants were determined recursively after a long 
computation exploiting the normalization requirements 
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FIG. 1. Pnum, which is obtained after 2 x 10 iterations of 
the map ( Jl5| ) , and Pthcor , which is the exact third order expan- 
sion in f3, are plotted for [3 = 10 (curves are indistinguishable). 
The inset compares /3(P num /Po — 1) to f3(P t i lcoI /Po — 1) (thick 
line). The Lyapunov exponent is A num ~ Athcor ~ 1.16293. 
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III. THE GOLDEN MEAN CASE: /3 = 1 
A. Generalities 

We now focus on the particular case j3 = 1 which ad- 
mits a non-perturbative treatment. An ingenious con- 
struction of the invariant distribution P(r) and the in- 
variant measure u([a, b]) — J a dr P(r) which involves a 
Stern-Brocot division of the real line has been proposed 
by Viswanath Jj7|. In this subsection we first recall that 
construction |j and the definition of the Stern-Brocot 
division (for details, see Ref. and also Ref. jl| which 
describes closely related Farey series). We then derive 
useful symmetry relations which will lead to new quanti- 
tative results concerning the invariant distribution P and 
the invariant measure v. 

The invariant distribution P{r) is symmetric, so we 
can limit ourselves to the half-line r > 0. The Stern- 
Brocot division is defined as follows: Start with the half- 
line Iq — and divide it into two intervals, 
and [r,h]', then divide every first generation interval and 
continue in this manner. Generally, the mediant of an 
interval [p/q,p' /q'] is obtained by "Farey addition": 



P 
<1 



p' p + p' 



(31) 



p p'~ 




p P + p' 


u 


p + p' p' 


q q'_ 




A q + q' _ 


q + q' q' 



Thus, the interval [p/q,p' /q'] is divided according to the 
rule 



(32) 



This generates 2™ intervals at the n th generation, each in- 
terval I of the (n — l) th generation producing a left and a 
right child (CI and 1ZI) in the n th generation, according 
to Eq. (p2|). For instance, [|, i] = TZCCIq, and gener- 
ally every Stern-Brocot interval I can be presented as 
/ = SIq, where S is the unique sequence of C's and IZ's. 
The length of S is denoted |«S|. Using this representation 
one can prove numerous properties of the Stern-Brocot 
division, e.g., the assertion that for every Stern-Brocot 
interval [|, |r], the identity p'q — pq' — 1 holds [§4|. 

The invariance condition, i.e. Eq. ( |l7| ) with 5 = 1, can 
be re-written in terms of the invariant measure to give 



< « 1 f ( 1 1 \ ( 1 1 

y,> 2 1 {-l + b'-l + aj \l + b'l + a 



Making use of the above equation, it is possible to ex- 
press the measure of the left and right children via the 
measure of the parent interval (?j : 



v{CSI ) 



(1 — r) v(SIq) if |<S| is even, 
r v(SIq) if |<S| is odd, 



(33) 



v{nSIo) -\(l-T)u(SI ) 



if \S\ is even, , , 
if 151 is odd, {6 > 



where r = (v5— 1)/2 is the golden ratio (which is also the 
inverse growth constant of the deterministic Fibonacci 
numbers, F n ~ t~"). 

We now turn to the derivation of symmetry relations 
which will be helpful later. To obtain the first one we take 
an arbitrary Stern-Brocot interval SIq = [|, ^7] and no- 
tice that the Stern-Brocot interval STIUIq = [|+2, ^+2] 
differs from SIq by a mere translation. Using Eqs.(|3 
(U) one finds r(l - r) v(SI ) = v(SKKI ), i.e., 



P „ v 
- + %—. 



q 



pv 



p p_ 

9' q' 



(35) 



with p = t(1 - r) = V5 - 2. 

Note that every interval with rational end points can 
be formed by joining appropriate Stern-Brocot intervals. 
Hence Eq. ( p5| ) holds for any rational interval [|, 2L]. 
Since rational numbers form a dense set on the line and 
v(r) = j/([0,r]) is a continuous (increasing) function, 
Eq. (|3^) applies to arbitrary intervals. Specializing to 
the interval [0,r], we arrive at the first relation 

v{[2,r + 2])=pv{r), (36) 

which can be also presented as 

P(r + 2)=pP(r) for r > 0. (37) 

The second symmetry relation, 

P(2 - r) = P[r), for < r < 2, (38) 

expresses a mirror symmetry with respect to 1 in the in- 
terval [0, 2]. To prove this result, we take a Stern-Brocot 
interval SCIq = [p/q,p'/q'] C [0,1] and notice that its 
symmetric image [2 —p'/q',2 — p/q] C [1,2] is also a 
Stern-Brocot interval. Specifically, the symmetric inter- 
val can be presented as SCIZI®, where S is the sequence 
obtained from S by exchanging the £'s and TVs. It turns 
out that the interval and its symmetric image have the 
same invariant measures 



v(SCIq) =v(scm 



(39) 



One then deduces Eq. ( p8| ) from Eq. ( J39| ) repeating the 
argument which has been employed in deducing Eq. (|37j) 
from Eq. ©. 

The proof of Eq. ([39]) can be accomplished by induction 
on the length. Suppose that Eq. ( ^9| ) holds for sequences 
S of length n. Assuming n even and taking the right 
child of SLIq, we use Eq. ( |34| ) to yield 



v (TISCIq) =(1-t)v (SCIq) , 

(cscnio) = (1 - r) v (scnio) 



(40) 
(41) 



and 



Since C = TZ and right-hand sides of above equations are 
equal, we obtain i/(1ZS£Iq) = ^(TZSCTZIq), thus prov- 
ing induction step for even n and the right child. The 
three other cases can be proved in a similar way. 
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B. Analytical results 

In this subsection, we will make use of the above 
symmetries to obtain quantitative results for the invari- 
ant distribution P and the invariant measure v. For 
instance, one can readily establish extreme behaviors. 
From Eqs. (|33|)-(|34|), one finds 



v{£ 2n I ) = v{n 2n Io)=p n v{Io) 



(42) 



for any (positive) integer n. Clearly, £ 2n Io = [0, 
and K 2n I = [2n,oo]. Additionally, z/(I ) = 1/2 which 
immediately follows from normalization. Thus Eq. ([42] ) 
reduces to 



From Eq. ( |43| ) we deduce the small r behavior, 

v{r) ~ e~ c /' with c = — — In p, 
and the large r behavior, 

\ v{r) ~ e~ cr . 



(43) 



(44) 



(45) 



Thus, the invariant measure v(r) sharply vanishes 
when r — > 0. Remarkably, similar plateaux exist around 
all rational r's. It is useful to demonstrate the exis- 
tence of the plateau at the proximity of a simple ra- 
tional point, say r — 1/2. Take the family of intervals 
£ 2n K£I = [|,f^]. The length of the n th interval 



is e„ = [2(4n + l)] -1 . Exploiting Eqs. fl33|)-(|34|), we can 
compute the measure of these intervals and then examine 
how the measure scales with length. We get 



1 2n- 
2' 4^~ 



1 



(l-rf 



c/4e 



(46) 



The continuity of the measure implies that the scaling 
holds for all small intervals, v(J^ + e) — ^(i) ~ e~ c / 4e . 

Similarly, one can work out the proximity of an arbi- 
trary rational number. One finds that the measure of 
the interval [p/q,p/q + e] (p and q are mutually prime 
integers) exhibits the following behavior 



exp 



q 2 e 



(47) 



with 



— \ m P as above. Similar asymptotics apply 
when we approach the rational point from the left. 

Turn now to the invariant distribution P(r). This func- 
tion obeys a striking number of intricate identities. [Note 
that it is a slight abuse of language to speak of P as a 
function: P is a distribution rather than a function.] In 
the following, we will be using Eq. @, with (5 = 1, so 
we write it down for easy reference 



2P(r- 



1 



P 



1 



(r + 2) 2 ~ V r + 2 , 
Note that Eq. (^8|) at r + 2n can be written as 

2P(r + 2n+l) = U n (r) + U n+1 (r), 
where we used the shorthand notation 
1 _ ( 1 



U n (r) 



P 



(r + 2n) 2 \r + 2n / 
Changing n to n + 1 in Eq. (|4^) gives 

2P{r + 2n + 3) = U n+1 {r) + U n+2 (r). 



(48) 



(49) 



(50) 



(51) 



Recalling that P(r + 2n + 3) - pP(r + 2n + 1) = [this 
is Eq. (37) at r + 2n], we reduce (|49|), ( |5l| ) to a recursion 
which involves only C/'s: 



U n+2 {r) + (1 - p)U n+1 {r) - pU n {r) = 0. 



(52) 



Since the variable r is mute, Eq. ( p2[ ) is merely a lin- 
ear recursion which is straightforwardly solved to find 
two independent solutions, (—1)™ and p n . Therefore, the 
general solution is 



1 



(r + 2n)'- 



P 



1 



2n 



A(r)(-l) n + B(r)p n . (53) 



In the n — > 00 limit, the left-hand side of Eq. (|53J) ap- 
proaches to zero. Thus, A(r) — 0. Another coefficient 
B(r) is found by specializing Eq. ([53|) to n — 0. Equa- 
tion (b3j) therefore reduces to 



1 



(r + 2n) : 



P 



1 



r + 2n 



(54) 



One can derive a few more useful formulae relating P(r) 
at different points. Performing the change of variable 
r — ► r _1 , one transforms Eq. (p4[) into 



P(r) 



(1 + 2nr) 2 



1 + 2nr 



(55) 



One can also take Eq. ( p4] ) at n = 1 and insert it into 
Eq. (Esl) . The outcome reads 



P(r + 1) = 



1 



2r 2 



"pt 1 - 



(56) 



Take now Eq. (|56|), change the variable r — ► r + 1, and 
use Eq. (]37f). This transforms Eq. (pq) into 



P(r) 



1+T 

1 + r 



P 



1 



1 



(57) 



We can take the same identity with (1 + r) _1 instead of 
r and insert it into the right-hand side of Eq. ( |57| ) thus 
obtaining another identity. Proceeding in this manner, 
we arrive at a series of identities 



G 



P(r) 



(l + r) 



2)n 



(F m +rF m _ l y 



P 



rF, 



m-2 



rF n 



(58) 



These identities apply for all integer m's as well as the 
Fibonacci numbers F m which are defined for all integer 
m, e.g., F_ 4 = 2, F_ 3 = -1, ^-2 - 1, F_i = 0, Fo = 1. 

Replacing r — ► 2 — r in identities ( |58| ) and using the 
symmetry relation (|3^) one obtains another infinite series 
of identities. The simplest such identity (corresponding 
to m = 1) reads 



P(r) = 



1+T 



1 



(59) 



The above identities together with the basic relation 
v([a, b]) = f dr P(r) can be used to re-derive all pre- 
vious results about the invariant measure, including the 
most interesting findings about the plateaux. 



interval (0, 1) onto the intervals which appear on the left- 
hand side. We then exploited Eq. ( p^ at n = 1, Eq. (p3S|), 
and Eq. (|57]), respectively, to obtain the second lines. 
Summing up Eqs. (|60])-(|62T) leads to identity 



(63) 



dxf(x)F(x)= f dx f(x)f[F(x)} 



where the linear operator T acts on F according to 



3 

f[F(x)]='%2a i F(h i (x)), 
i=i 

3 

Oil = Pi a 2 = a 3 = (1 + t)~ 2 , y^ctt = 1, 

i=l 



/ii(a;) 



.7" 



1 + 2a; 



, h 3 (x) = 



3 — x 



1 + x 



(64) 

(65) 
(66) 



C. Generating the invariant measure 

Because of symmetries, it is sufficient to study v{r) on 
the interval [0,1]. To distinguish the restricted version 
from the general case we change the notation: r — > x, 
/, f ^r-N. The prefactor (1 - p)/4 guaran- 
tees that the invariant measure N(x) = dy f(y) obeys 
N(l) = 1. To see this, we take the normalization condi- 
tion, 1 = J_ drP(r), and massage it a bit: 



r°° 2 f 2 4 f 1 

2/ drP(r) = / drP{r) = / dr P 

Jo l-PJo 1 - P Jo 



(r). 



The first above identity is obtained by using Eq. j37| ) 
and summing up the geometric series while the last is an 
obvious consequence of Eq. d38|). Thus we indeed obtain 
N(l) = 1 if we choose f(x) = P(x). 

Now let us integrate f(x)F(x) with any F(x) over in- 
tervals (0, 1/3), (1/3, 1/2), (1/2, 1). We find 



dx 



1/3 rl 

dxfF = Jo {l + 2xf J \l + 2x 



f 



F 



p / dxf(x)F , 

' Jo ^ ' V I -2.r 

I- 1 



X 



1/2 r dx 1 1 , 



1/3 



3 — x 



1 + 2x, 
(60) 

1 



3 — x 



(l + r)" 2 / dxf(x)F 







3 — x 



, (61) 



dx 



f 



1 



1/2 Jo (1 + ^) 2 VI + x 

= (l + r)- 2 / dxf(x)F 



F 



1 + x 



1 + x 



(62) 



In deriving the first lines in above formulae we have used 
the mappings :r(l + 2a:)~ 1 , (3 — a;) -1 , (1 + a;) -1 of the unit 



These relations show that the distribution / can be gener- 
ated by the following simple recursion process. We take 
an arbitrary point x € (0,1] and assign a unit weight 
w(x) = 1. We then define three new points and associ- 
ated weights according to the rule 

Xi = hi(x), w(xi) = atiw(x), i = l,2,3. (67) 

Iterating n times, we get 3™ different XiS. The weights 
w(xi) of these points add up to 1 and w(x) converge (in 
the weak sense) to f(x). In Fig. 2, we plot the invariant 
measure obtained after 4 and 10 iterations. 



2 



W(3/5)=6x-3- 



N(4/5)=8t-4- 

-N{1/2)=t 
-W(2/5)=4x-2 



-N(3/4)=4-5t 
-N(2/3)=2t-1 



— N(1/3)=3-2t 

N(1/4)=2-3x 

-N( 1/5) =5-8t 



FIG. 2. The invariant measure N(x) — J dy f(y) obtained 
after 4 iterations (dotted line involving only 3 4 = 81 different 
Xi's) and 10 iterations (3 10 = 59049 different xt's). 

The plateaux described by Eq. (|4^) are clearly vis- 
ible in Fig. 2. The smaller the denominator q, the 
wider is the associated plateau as predicted by Eq. (|47|) . 
In a pseudo-gap of / (associated to a plateau for N), 
N{p/q) = {n p / q T}, where {.} denotes the positive frac- 
tional part and n p / q is a possibly negative integer. The 
values of N(p/q) for a few small denominator fractions 
are presented in Fig. 2; they have been calculated by us- 
ing Eqs. (|63|)-(|66|). Thus, a kind of "pseudo-gap labeling 
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theorem" p5[ applies here: for tight-binding Hamilto- 
nians on ID quasiperiodic chains associated to an irra- 
tional number a (see Ref. p5| for a general formulation), 
this theorem states that the integrated density of states 
(IDOS) displays plateaux (corresponding to the energy 
true gaps). In these plateaux, the IDOS takes values pre- 
cisely of the form {na}. In one dimension, the IDOS is 
intimately related to the invariant measure v [0. There- 
fore, it is not surprising that v exhibits a similar behavior. 

We now show how to calculate the Lyapunov exponent 
from the restricted invariant distribution f(x). Take the 
basic formula, A = J dr P(r) lnr, and perform the same 
massage as we did for the normalization condition in the 
beginning of this subsection. We obtain 

A = / dxf(x)g(x) 



1-P 



1 -x 



2 ^ - — y V 2fc + 1 

k=0 x 



= (1 - p)J2 P k ln ( 2fc + 1 ) = 0.29320491137 . . 

k=Q 

We then numerically compute A via 

3" 



A 



= lim w(xi)g(xi) + p. 

n — >oo * — • 



For instance, starting from x = 1 — r, one reproduces 10 
digits of A after only 8 iterations. 



IV. GAUSSIAN RANDOM SEQUENCES 



1 



(3z v 



Rn-l 

The invariant distribution P(R) satisfies 



(70) 



P(R) = J dz G(z) J dR' P(R')S ^R-l-^j, (71) 



which can be re-written as 



P(R + 1) = — L= / dR' \R'\ P(R') e 



(72) 



Note two general features of the invariant distribution 
which directly follow from Eq. (|72|). One is a power-law 
large R asymptotics, 



P(R) 



2 (3P(0) 



for R > /3. 



(73) 



Another unexpected feature is a weak logarithmic diver- 
gence at R = 1: 



for \R-1\ < 13. (74) 



Equation (||) follows from Eqs. @-(|73|). 

We now turn to a perturbative analysis. While for 
the Gaussian random recurrence, there appears to be 
no threshold separating weak and strong disorder, dif- 
ferent approaches should be implemented when [3 — > 
and (3 — > oo, respectively. In the former region, the regu- 
lar perturbation theory applies while in the latter region 
one needs a singular perturbation theory. 



Invariant measures exhibited by the Fibonacci random 
sequences appear fractal for all strengths of disorder. 
This is apparently caused by the discreteness of the dis- 
order. Additionally, there is a transition between weak 
disorder, manifested by an invariant measure whose sup- 
port is bounded, and strong disorder characterized by an 
invariant measure with unbounded support. This fea- 
ture is seemingly caused by the boundness of the disor- 
der. Hence, random sequences with unbounded continu- 
ous disorder distributions provide a natural counterpart 
to random sequences with binary disorder. As a specific 
example, we consider the Gaussian random recurrence, 



x n -\.\ x n -\- [3z n x n ~\ 7 



(68) 



where z„'s are independent identically distributed ran- 
dom variables with Gaussian probability density 



G{z) 



1 



(69) 



To analyze the Gaussian random recurrence, we again 
use the Ricatti variable R n — x n+ i/x n which reduces 
Eq. (|68|) to the random map 



A. Weak disorder 

Equation (|7(]) shows that for weak disorder (/3 — ► 0), 
the magnitude of R — 1 is comparable with (3. This sug- 
gests to rescale the variable R and the invariant distri- 
bution P(R) according to 

R = l + 0r, P(R) = p- l Q(r). (75) 

The normalization condition J dR P(R) = 1 now reads 

J drQ(r) = 1. (76) 

The governing equation (|72] ) becomes 



Q(r) 
G(r) 



dr'\l + (3r'\ Q(r') e -- 2 (^'+/3V 2 /2)_ (77) 



In the small /3 limit, we simplify Eq. (|77]) to 



Q(r) 
G(r) 



dr' (1 + f3r') Q(r') e -^^'+^^l^. (78) 
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The error caused by the above simplification is of order 
Q{—l/(3). We will see that it vanishes as exp(— 1/2/3 2 ) 
and therefore it can be ignored in perturbative analysis. 

We are seeking a perturbative solution. The symmetry 
(3 «-* — (3 suggests an expansion in (i 2 rather than (3: 



(79) 



n=0 



We must also expand the exponent on the right-hand side 
of Eq. ( fr"7| ) . This exponent is the generating function of 
the Hermite polynomials H n (r), 



cxp[-r 2 (a; + .T 2 /2)] = 



where H n {r) are defined via EE 



(—xr) T 



H n {r), (80) 



(81) 



Thus, the Hermite polynomials make their appearance in 
the problem with non-Hermitian Hamiltonian. 
Inserting Eqs. (f79|)-([30"|) into Eq. (f78|), we obtain 



Qn(r) 
G{r) 



£ 

fc=i 



„2fc-l 



(2k)\ 



H. 



2k+l\ 



(r) / dr'(r') 2k Q n _ k (r'). (82) 



In deriving Eq. (|82|), we also used the recursion relation 
for the Hermite polynomials f23|]: 



H n+ i(r) = rH n (r) - nH n ^x(r). 
Solving Eq. recursively yields 



(83) 



Qo(r)=G(r), 
Qi(r) 1 



G(r) 
Q 2 (r) 
G(r) 
Q 3 (r) 
G(r) 

Qi{r) _ 7857 
G(r) ~ ~32~ 
7 

+ 16 



rH 3 (r), 



^ 7 :rH 3 (r) + -r i H 5 (r), 



12rH 3 (r) + ^ r 3 H 5 (r) + i r 5 H 7 (r) 
1 35 

rH 3 (r) + ^r 3 H 5 (r) 



384 



r 7 H 9 (r), 



Q 5 (r) 1362843 48675 3 

-^y = o Cg ^(r) + -tt^ r #s(r) 



256 



128 



-— r 5 H 7 (r) + — r 7 H 9 (r) + — 
16 n ' 32 9W 3840 



'^ii(r), 



etc. We will use these results to compute the Lyapunov 
exponent. The basic formula (ph now reads 



A = 



drQ{r) ln|l + /3r|. 



(84) 



Expanding the logarithm and the invariant distribution, 
Eq. (f79l), and recalling that Q(r) — Q{—r), we get 



00 ft -, /.oo 

A =-E^"EW drr 2fe g„_ fc (r). 

n=l fc=l ^° 



(85) 



Inserting above expressions for Q n (n = 0, . . . , 5) into 
Eq. (p5|), we obtain the weak disorder expansion: 



22 f3 e 



13197 



32 



6 s 



2374335 
256 



/5 



10 



118392093 



512 



/3 12 +0(/3 14 ). (86) 



Interestingly, neither the R~ 2 asymptotics ((73|) nor the 
logarithmic singularity ( [74] ) appear in the weak disorder 
expansion. Both these behaviors are non-perturbative. 
For instance, P(0) which appears on the right-hand sides 
of Eqs. d73])- ( |74|) scales as exp(— 1/2/3 2 ), i.e., the behav- 



iors Eqs. (73)-(|74|) are beyond the scope of perturbation 
techniques. Note also that the observation of the loga- 
rithmic singularity (Q) requires probing a prohibitively 
tiny region r ~ exp[— exp(— 1/2/3 2 )]. 

The occurrence of non-perturbative corrections (of or- 
der exp(— 1/2/3 2 )) suggests that the radius of convergence 
of the series ( |S6| ) is equal to zero. This is (non-rigorously) 
confirmed by the following approximate analysis of A(/3). 
First, we write R n « 1 + (3z n , as the distribution of R n -i 
is strongly peaked at R = 1 for small (3. Within this 
approximation, and using the fact that the Gaussian dis- 
tribution is an even function, we obtain 



dzG{z)ln\l- (3 2 z 2 \dz, 



(87) 



which reproduces exactly the first term of Eq. (|8^) . Now, 
expanding Eq. (^) in powers of 1 , is likely to lead to 
the correct qualitative behavior for the full general ex- 
pansion. The generic term, 



In r°° 



2n 



dz G(z) z 



In 



-rU + - 



grows faster than any exponential, ensuring that the ra- 
dius of convergence is indeed zero. Of course, we cannot 
exclude that for the actual expansion, subtle cancella- 
tions lead to a finite radius of convergence. However, the 
occurrence of non-perturbative corrections, and concrete 
ingredients of the argument presented above (mainly, the 
unboundness of the distribution G(z) and the fact that 
the series ln(l+a;) as a finite radius of convergence) which 
seem to persist in the general case, favor a zero radius of 
convergence and the asymptotic character of the series 
(|S6|). This is very different from the case of random Fi- 
bonacci sequences where the weak disorder expansion has 
a finite radius of convergence and there were no trace of 
any non-perturbative contribution. 
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B. Strong disorder 

For P — > oo, we again use the properly normalized Ri- 
catti variable y n = x n +i/x n \f]3. Equation ( |68| ) reduces 
to the random map 



Vn = — + S, 6 = p 

Vn-1 

The invariant distribution satisfies 
1 



-1/2 



(88) 



P(y + 6) 



'2tt J- 



d^rjlPfaJexH-O-^ (89) 



and the Lyapunov exponent is given by Eq. (|16| ) as in the 
Fibonacci case. Equation ( {39| ) suggests to seek a pertur- 
bative solution. In the zeroth order, one might set 5 = 
in Eq. (|S9"|). The corresponding invariant distribution 
Pq (y) is an even function of y which satisfies 



Po(y) 



dr] r] Pq(t]) exp 



2 2 



Paradoxically, a (formal) solution to this equation, 



(90) 



(91) 



does not obey the normalization requirement. This indi- 
cates that the naive perturbation approach does not work 
and one must develop a singular perturbation theory. 
One still anticipates that Po(y) is given by Eq. (|9l]) apart 
from the small and large scales, \y\ ~ S and \y\ ~ 
which are implied by the random map (p^). Treating 
these scales as cutoffs allows us to normalize the solu- 
tion ( pi] ) and to estimate the amplitude A » (2 In/?) -1 . 
One can establish the existence of the cutoffs more rigor- 
ously. Using Eqs. @, @ one estimates P(0) - A5 -1 
in agreement with the existence of the small scale cut- 
off y ~ S. The large scale cutoff already follows from 
Eq. @ which now reads 



P(y) 



2 P(0) 



7T 



for y ^ S 



(92) 



Note also that the deficiency of the naive perturba- 
tion approach is clear from the respective random map, 
y n = z n /y n _i. Indeed, iterating the above map and tak- 
ing the logarithm gives In y n = -0" ~ fe m z k- The cen- 
tral limit theorem now asserts that the scale of the lim- 
iting the distribution grows indefinitely with n, namely 
In y n ~ y/n. Thus for 5 = already the basic assumption 
that y n approaches to a limiting distribution which does 
not depend on n is incorrect. 

We now present a computation of the zeroth order con- 
tribution to the Lyapunov exponent which does not re- 
quire the knowledge of Pq. Denote by A the zeroth order 
contribution to the Lyapunov exponent. We have 



/>oo 

A = 2 / dyP (y) Iny 
Jo 

= 2y-y drjr)P (ri) J Iny exp 



2 2 

V V 



dr)P (r)) 



*( 2 ) +ln2-21nr/ 



!• I - I In 2 



A. 



In the second line we used Eq. (|90|) ; this step is not really 
rigorous though we think the final result is correct. The 
variable t which appears in the third line has been defined 
via t = y 2 f] 2 /2; in the fourth line we used the digamma 
(psi) function, V(x) = T'(x)/T(x) f23|; in the last line we 
used the normalization requirement and the definition of 
A. The above equation yields A which can be simplified 
further by using the identity $(1/2) = —7 — 2 In 2, where 
7 is the Euler constant. Finally, 



A = - 



7' 



In 2 



= -0.317590711365. 



(93) 



Our numerical results suggest that the strong disorder 
expansion involves powers of (ln/3) -1 rather than (3~ l : 



1 



A(/3) = -ln/3 + A + ]TMln/3)- 



(94) 



fe=i 



Of course, it is hardly possible to probe higher order log- 
arithmic terms numerically. However, plotting A — i ln/3 
versus (ln/3) -1 gives a fairly straight line for (3 > 10 4 , 
with the slope b\ w 0.557, and a perfect fit to the above 
functional form, keeping a quadratic term in (ln/3) -1 , 
with b 2 ~ -0.52 (see Fig. 3). 




FIG. 3. We plot A(/3) for Gaussian random sequences. The 
upper insert compares -(A(/3) + f3 2 /2)/f3 4 = 9/4 + 22/3 2 + ... 
obtained from the 10th and 12th order expansions in Eq. (86) 
to the result of numerical simulations. The lower insert is a 
quadratic fit in (ln/3) -1 (see Eq. (94)) of A(/3) — \ In (3, for large 
j3, leading to an extrapolated A ~ —0.317, with 61 » 0.557 
and 62 « -0.52. 
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Extrapolating the quadratic fit to (ln/3)^ 1 = yields 
A w —0.317, in good agreement with the theoretical pre- 
diction Contrary to the Fibonacci case, the curve 
of A(/3) appears perfectly smooth, and is certainly not 
fractal. 

Finally, we note that asymptotic methods |23] should 
in principle allow to perform the strong disorder expan- 
sion more systematically. 



V. DISCUSSION 

The rich behavior exhibited by random Fibonacci num- 
bers suggests avenues for further investigation. For in- 
stance, how to reconcile perturbative results in the large 
(3 limit with non-perturbative results for (3 = 11 This 
question is important as there appears to be just a sin- 
gle threshold (3 = 1/4 and therefore f3 — 1 lies within 
the large (3 domain. This suggests qualitative similar 
behaviors which is not the case. The major difference 
between (3 = 1 and (3 — > oo cases is manifested in ex- 
treme behaviors of the invariant measure. In the former 
case, it exhibits exponential asymptotics, Eqs. (Q)-([l5]), 
while the latter is characterized by power-law asymp- 
totics, v(r) r for r — > and ^ — v{r) ~ for 
r — > oo. More generally, our perturbative results are 
infinitely smooth, in a gross disagreement with the be- 
havior for (3=1. Another (related) set of questions con- 
cerns the curve X((3): Is it a fractal? Does it become 
genuinely smooth at least for sufficiently large (3 ? 

One could ask for a more complete characterization of 
the growth (decay) of the random Fibonacci numbers. A 
natural conjecture is x n ~ e Xn n w . For the neutrally sta- 
ble recurrence, x n +\ = x n ± (3 s x n -\ with (3 a w 0.70258 
chosen so that X((3 S ) = 0, the above conjecture would 
imply the power-law behavior x n ~ n" s . 

Finally, it would be very interesting to analyze random 
Fibonacci numbers for the critical strength of disorder, 
(3 = 1/4. This strength of disorder appears a bit more in- 
teresting than (3=1: It is hard to see what distinguishes 
(3=1 from say (3 = 1.23456, while the case of (3 = 1/4 is 
certainly special as it separates the regions of weak and 
strong disorder. 

We are grateful to H. Castillo, G. Oshanin, S. Redner, 
and P. Sharma for discussions. The work of PLK was 
supported by NSF Grant DMR9978902. 
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